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Abstract. - We show that a new glassy phase can emerge in presence of strong magnetic frus- 
tration and quantum fluctuations. It is a Valence Bond Glass. We study its properties solving 
the Hubbard-Heisenberg model on a Bethe lattice within the large N limit introduced by Affleck 
and Marston. We work out the phase diagram that contains Fermi liquid, dimer and valence 
bond glass phases. This new glassy phase has no electronic or spin gap (although a pseudo-gap 
is observed), it is characterized by long-range critical valence bond correlations and is not related 
to any magnetic ordering. As a consequence it is quite different from both valence bond crystals 
and spin glasses. 
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The interplay of strong quantum fluctuations and ge- 
ometrically frustrated magnetic interactions can give rise 
to new low temperature phases. As noticed by Ander- 
son [1] a way to minimize the effect of frustration and 
obtain a low energy state is coupling the electrons in va- 
lence bonds. A very good variational wave function that is 
generically in competition with the antiferromagnetic (or 
more general magnetic) state can be obtained by form- 
ing a superposition of short range valence bonds that are 
arranged as dimers on the lattice. If no lattice symme- 
try is broken this corresponds to the (so called) resonat- 
ing valence bond liquid (RVBL). In the last decades, this 
state has received a lot of attention in connection with the 
unusual physical behavior of the normal phase of under- 
doped high T c superconductors [2]. Indeed Anderson [3] 
proposed that the holes created by doping the antiferro- 
magnetic insulator (of the high T c 's phase diagram) can 
gain substantial kinetic energy in the RVBL state and not 
in an antiferromagnetic background. As a consequence, 
doping favors the RVBL state which could then become 
the thermodynamic stable phase and be responsible for 
the unusual behavior of underdoped samples. Concomi- 
tantly, resonating valence bond ground states have been 
the focus of an intense activity [4] in the context of frus- 
trated magnets. RVBL or spin liquids have been found 
for several models [4]. These states can undergo quan- 
tum phase transitions where lattice symmetries are spon- 
taneously broken. This gives rise to valence bond crystals 
(VBC). Different models are known to lead to this type of 
ground states [4] characterized by long range dimer-dimer 



correlations. The situation in experiments is complicated 
by unavoidable magneto-elastic couplings: making the dif- 
ference between induced and spontaneous dimcrization is 
a difficult task. A first experimental example of sponta- 
neously broken states has been apparently found in [5], 
The aim of this work is to study a new kind of valence bond 
state: the valence bond glass (VBG). Similarly to VBC the 
arrangement of the dimers (or valence bonds) breaks the 
lattice symmetry. However, contrary to VBC, this corre- 
sponds to an amorphous dimerization and not crystalline 
one. Although VBG arc analogous to spin glasses [6] they 
are physically quite different. In particular the spins do 
not freeze in a disordered profile. We expect that the VBG 
phase can arise in presence of strong magnetic frustration 
as one of the competing ground states. The addition of 
(little) quenched disorder will favor this phase. Depending 
on the system, the low temperature phase could be either 
a VBG or a spin glass. Actually, the spin glass phase is 
conjectured to exist even in absence of disorder on some 
frustrated lattices [7] (see however [8,9]). 
In the following we shall investigate the properties of 
the valence bond glass phase focusing on the Hubbard- 
Hcisenbcrg model within the large N approximation intro- 
duced by Affleck and Marston [10]. The underlying lattice 
we shall focus on is a random regular graph with connec- 
tivity z 1 . The reason for this choice is twofold. First, 
this type of graphs are on any finite lengthscale as Bethe 



It is a graph taken at random within the set of graphs whose Af 
sites are all connected to z randomly choosen neighbors. 
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lattices or Cayley trees. This, as it is well known for clas- 
sical systems [11], introduces useful simplification in the 
analysis of the model. The main reason is, however, that 
topological frustration and quenched disorder are intro- 
duced by very long loops (of the order log TV where N 
is the number of sites) in random regular graphs. These 
loops disfavor crystalline states and let emerge easily the 
glassy phases [12,13]. We consider the SU(N) version of 
the Hubbard-Hciscnbcrg model introduced in [10]: 
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where Cj iQ denotes the destruction operator of an electron 
of spin index a (a — 1, . . . , N with N even) on the site 
i. The sum is restricted on nearest neighbor sites on 
the lattice. The first two terms correspond to the SU(N) 
Hubbard model, while the last term accounts for the near- 
est neighbor antiferromagnetic interaction (J > 0) 2 . We 
shall focus on the N — > oo limit and consider only the half- 
filling case, where rn/(N/2) = E a 4, a *,*/( N / 2 ) = 1 for 
all sites. Using that Si&j equals - Y^ a ,§ ^a^i,^ up 
to constant terms in the large N limit [10, 14], the Hamil- 
tonian can be rewritten in a SU(N) manifestly invariant 
form. At half-filling it reads: 
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Note that all terms constant in the large N limit have 
been neglected. Here and henceforth the summation over 
the SU(N) indices will be skipped for simplicity. The 
partition function of the system at finite temperature can 
be written as a path integral 
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where (3 is the inverse temperature, and the (imaginary 
time) Lagrangianis £(c, c^) = T~C+J2i c la (d/dr)cj i0 . The 
functional integral is of course non trivial, due to the pres- 
ence of the non linear interaction. However, one can per- 
form a Hubbard-Stratonovich transformation which allows 
to rewrite the Lagrangian quadratically in the fermions, at 
the expense of introducing a new (complex) bosonic field, 
Xij, on each edge of the lattice [10]: 

£(c,c f ,x) = 5Z{-dXtf| 2 - (* + Xy) c !,<* c j> + h - c 
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The equation of motion of the auxiliary bosonic field reads: 



J 



(Xi j (r))= ^<4aWd,a(T)). 
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Xij is the valence bond field and gives an extra contribu- 
tion to the electron hopping amplitude between the sites 
i and j. The number of valence bonds on link (ij) is given 
by N\xij\ 2 /J up to subleading terms [10]. 
The advantage of this representation is that the integral 
over the fermionic degrees of freedom is now Gaussian. 
Therefore, they can be integrated out, leading to an effec- 
tive action which depends only on the bosonic variables: 



exp[-S e ff(x)] = JvcV^exp - J AtL{c,c\x) 



(6) 



The effective action thus reads: 



rP i 

S eff =N dr -r \XH I 2 - ^Tr log M, (7) 
Jo dd) J 

where the matrix M is given by M = [(d/dr)I - tC - x], C 
being the connectivity matrix of the lattice, i.e., C,j = 1 if 
i and j are nearest neighbors on the lattice and zero other- 
wise, x h as an analogous definition except that Xij = Xij 
if i and j are nearest neighbors. 

So far, these transformations are exact and do not depend 
on the particular choice of the lattice. In the TV — ► oo 
limit the saddle point integration over the bosonic vari- 
ables, Xij i becomes exact and we can compute the free en- 
ergy of the system by seeking the lowest minimum of the 
effective action 3 . Assuming that at the saddle point the 
valence bond operators are time-independent, the prob- 
lem reduces to finding the minima of the "classical" free 
energy (3F(x) = S eff /N (N being the number of SU(N) 
indices), 
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We denote by A the eigenvalues of the one-particle Hamil- 
tonian 



W i = -E [(t+xij)4 Cj +h, 
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Note that the (complex) bosonic variables Xij can have any 
arbitrary spatial dependence and that there is no need to 
introduce the chemical potential since it is expected, and 
found, to be zero at half filling 4 . For simplicity we will 



(4) 



2 As discussed in [10], the antiferromagnetic interaction is not 
generated in perturbation theory at N = oo, so it has to be added 
in the original Hamiltonian. 



3 If we had decoupled the U term in eq. 1, as done for the J term, 
by introducing a field 0; then we would have found saddle point 
equations leading, at half filling, to the solution <j>i = [10]. That is 
the reason why we dropped this term from the beginning. 

4 Although random regular graphs arc not bipartite, they behave 
in a similar way. In particular, for all phases, we find electronic 
densities of state that are symmetric around zero. Thus, the chemical 
potential is zero at half filling. 
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set J = 1 in the following, bearing in mind that all energy 
scales are measured in units of J. 

The saddle point equations consist simply in Eq. (5) where 
the average on the RHS is performed using the Hamilto- 
nian "ri\. Obtaining an analytical solution for a given par- 
ticular lattice is, in general, a hard task. However, in some 
special cases, the problem can be simplified. In particular 
by considering periodic solutions one reduces the indepen- 
dent degrees of freedom to a finite number (4 in the case 
studied by Affleck and Marston [10]). Our aim is to find 
whether there are amorphous or chaotic solutions. Thus, 
in our case, obtaining a full analytical solution seems ex- 
tremely difficult. 

On infinite random graphs the Bcthc-Pcicrls approxima- 
tion is exact [12]: since the average length of the loops 
is infinite, it is possible to write down self-consistent it- 
eration equations for local "cavity" Green's functions, 
(or "Weiss functions"), Gi, defined on each site of the 
graph [15]. In particular, for any given configuration of 
the valence bonds, {xij}> it is straightforward to show 
that the following recursion relations must hold: 



(10) 



where v n = [2n + l)7r//3 are the fcrmionic Mat- 
subara frequencies. The Green's function, Gi(v n ) = 
—/3(ci a (vn)Cia(vn)}, can be calculated on each site as 
a function of the Gi on the neighboring sites, by using 
Eq. (10), where the sum is extended over all the z neigh- 
bors. For any given finite graph, and for any given profile 
of the bosonic field, Eqs. (10) provide a set of solvable 
equations for the cavity propagators. Furthermore, by 
enforcing the equation of motion for the valence bonds, 
Eq. (5), one finds that, on each link of the graph, the 
bosonic operators must verify: 



Xij 
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The last equation is non-local, and is reminiscent of the 
TAP equations derived in the context of spin glasses [17]. 
For infinite systems Eqs. (10) and (11) allow to treat the 
liquid and the dimer phase (see below) in a very natural 
way. The analysis in the glass phase is much more involved 
and complicated. See [12] for the method used in classical 
cases 5 and [16] for its extension to quantum cases. As a 
consequence we will use the previous approach to study 
simple (non disordered) phases and the transition lines. 
In order to study the glassy phase we interpret the free 
energy, Eq. (8), as the Hamiltonian of a classical system 
of complex variables. Hence, the problem of finding the 
minima of the free energy is reduced to finding classical 



5 The cavity method that would be needed to analyze the glassy 
phase is substantially more difficult than the one developed for spin 
glasses on Bcthe lattices. The reason is that the valence bond inter- 
action is on all scales and not only between nearest neighbors. 




Fig. 1: Phase diagram of the Hubbard-Heisenberg SU(N) 
model at half filling on the random regular graph (z = 3). We 
show the relative positions of the uniform phase (U), the dimer 
phase (D), and the valence bond glass (VBG). At t a (T) the 
uniform phase becomes unstable, the valence bond non-linear 
susceptibility diverges (see Fig. 3), and a continuous transition 
from the Fermi liquid to the VBG takes place. At t c (T) the free 
energies of the dimer phase and that of the VBG coincide and 
a first-order transition occurs. The dashed line corresponds to 
the spinodal of the dimer phase. The probability distributions 
of the valence bonds in the different phases are reproduced 
schematically in the insets. 



ground states. To solve the latter problem we use Monte 
Carlo annealing simulations. Basically, we introduce an 
auxiliary temperature T aux and, at each step, we attempt 
to change one Xij a t random according to the Boltzmann 
weight e _F ( x )/ Taux . The move is accepted with probability 
p = min {1, exp[— AF/T aux ]}. The auxiliary temperature 
is finally decreased at constant rate down to zero temper- 
ature. Details on the numerical procedure are discussed 
in the Appendix. 

By employing both the analytical and the numerical ap- 
proaches described above, we have derived the phase di- 
agram of the SU(N) Hubbard-Heisenberg model on the 
random regular graph with connectivity 2 = 3, see Fig. 1. 
Uniform phase — At high enough temperature and hop- 
ping amplitude the system is in a uniform phase, where 
the bond operators are real and equal on each link of the 
graph, Xij — X- F° r a given value of x, the electronic 
density of states can be computed easily since the density 
of states of the connectivity matrix is known [18], see the 
inset of Fig. 2. The uniform phase is translational invari- 
ant and gapless. It is clearly a Fermi liquid. 
For each value of T and £, x{T, t) in the uniform phase can 
be computed within the Bethe approximation, by enforc- 
ing translational invariance into Eqs. (10) and (11) (i.e., 
Gi = G and Xij = x)j which reduce to a simple algebraic 



Marco Tarzia and Giulio Biroli 



equation: 

* = E 



(t + x)/P 



it + (2,-2)1* ■ 



v -f + {z-l)\t- 



(12) 

One can then check the stability of the liquid solution 
with respect to any other solution of the bosonic field. 
This amounts in studying the (lowest) eigenvalues of the 
Hessian of F(x)- Using the base where the one-particle 
Hamiltonian, Eq. (9), is diagonal, and Fourier transform- 
ing with respect to the imaginary time, one gets 



1 _ e /3(A+Y) 



A, A' 



A + A' 



(13) 
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+ ef 3x + eP x ' + eA x+x ') ujI + (A + A') 2 



where v\ is the z-th component of the eigenvector asso- 
ciated with the eigenvalue A, and uj n = 2m: / (3 are the 
bosonic Matsubara frequencies. The first instability of 
the uniform solution is expected to correspond to a long 
wave-length modulation and should thus occur at uj n= o 
first. In order to analyse it, we generate random regular 
graphs of size Af and compute A, v\. Then using Eq. (13), 
we find that the smallest eigenvalue of the Hessian ma- 
trix at zero frequency becomes negative as either T or 
t are decreased down to t s (T,J\T). We then extrapolate 
the value of t s (T, Af) (averaged over several realisations of 
the graph) in the Af — » oo limit by increasing Af from 64 
to 1024. The curve t s (T) in the thermodynamic limit is 
shown in Fig. 1. In particular, at T — the liquid solution 
becomes unstable at t s ~ 0.29. 

Dimer phase — At low enough temperature and hopping 
amplitude a dimer phase (or Peierls phase) [10] is found to 
minimize the system free energy. In this phase the valence 
bonds can assume only two possible values, Xi on N/2 
links and xi on the others Af(z — l)/2, with \\\ \ > |x2 1 , in 
such a way that each site has exactly one link where the 
bosonic operator equals \i an d z — 1 links where it equals 
Xi- As the random regular graph is dimerizable [19], the 
analysis of Ref. [14] guarantees that a dimer phase (with 
Xi = 0) is the actual ground state of the pure antiferro- 
magnetic system (t = 0). 

At any given temperature and hopping amplitude, xi an d 
Xi can be determined analytically within the Bethe ap- 
proximation. More precisely, one allows the cavity Green's 
functions and the valence bonds to assume only two possi- 
ble values, respectively Q\ and Q2, and xi an d Xi- Taking 
into account the structure of the dimerized configurations, 
one can obtain a closed set of equations, which can be eas- 
ily solved: 



|t±X2 



Q a (y n ) 

Xa(b) 



1 t + Xa(b) 



\t±Xi 



if a = 1 
if a = 2 



13 n [06(a) ("n)] - I* + Xa{b)\' 



(14) 




Fig. 2: Main frame: Overlap probability distribution, P(q), 
at zero temperature and t = 0.23 in the VBG. The data are 
averaged over 16 different realizations of the graph, with M — 
256. The delta function in q = corresponds to the fraction 
of replicas which end up in the same state, and it is expected 
to disappear in the thermodynamic limit (e.g., for a system 
of AT = 128 sites the delta peak in zero is approximately 1.5 
bigger than that for Af = 256). Inset: Electron spectrum, 
jo(A), at zero temperature in the different phases: Fermi liquid 
at t — 0.34 (corresponding to the point marked by x in the 
phase diagram of Fig. 1, dotted line), Valence Bond Glass at 
t = 0.23 (point marked by • in Fig. 1, continuous line) and 
Dimer phase at t — 0.16 (point marked by ■ in Fig. 1, dashed 
line). The electron spectrum has been computed analytically 
in the uniform and in the dimer phase, and numerically in the 
VBG phase. 

In the dimer phase, both xi an d X2 turn out to be real 
(but at t = 0, where the system has a local gauge sym- 
metry, da — > Ci a e l8i and c\ a — > c\ a e~ iei ). The electron 
spectrum in the dimer phase can be found similarly by 
computing the resolvent of the matrix t<C + x in the dimer- 
ized state. The (electronic) density of state has gap, see 
inset of Fig. 2. This also induces a gap in the spin exci- 
tations 6 . Using the above results, the free energy of the 
dimer phase can be determined exactly for each value of 
T and t. At small enough temperature and hopping am- 
plitude the dimer phase corresponds to the absolute min- 
imum of the free energy. For larger values of t (or T) the 
dimer phase reaches the spinodal line, where the gap closes 
and the smallest eigenvalue of the free energy Hessian ma- 
trix vanishes (dashed line in Fig. 1). At zero temperature 
this happens at t ~ 0.218. Note that this zero temper- 
ature spinodal point lies below the corresponding one of 
the liquid which is the stable phase at high t. As a conse- 
quence, there is necessarily an intermediate phase. As we 
shall show in the following this is the Valence Bond Glass. 

Valence Bond Glass — In order to study and prove the 
existence of the Valence Bond Glass phase we use Monte 
Carlo annealing simulations for the reasons explained pre- 



6 The spin Green function can be obtained quite easily from the 
electron Green function in the large N limit [10]. 
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viously. First, we check that our numerical procedure gives 
back the uniform (dimer) phase at high (low) enough tem- 
perature and hopping amplitude. In the intermediate re- 
gion where both phases are unstable (e.g., at zero tem- 
perature for 0.218 < t < 0.29) we find that amorphous 
configurations of Xij correspond to the actual minima of 
the free energy. This is a glassy phase, which we call va- 
lence bond glass. This is not a spin glass since the average 
value of the spin is zero on each site of the lattice, (Si) = 0, 
as the SU(N) symmetry is unbroken. 
The valence bonds, Xiji are rea l valued and their dis- 
ordered profile is described by a nontrivial distribution, 
P (x), as schematically depicted in the inset of Fig. 1. The 
electron spectrum is gapless in the VBG, although it ex- 
hibits a pscudo gap, as shown in the inset of Fig. 2, which 
becomes deeper and deeper as either the temperature or 
the hopping amplitude are decreased. 
Interestingly enough, similarly to spin glasses [6], on any 
given graph different annealing procedures may lead to dif- 
ferent configurations with the same free energy. One can 
measure the distributions of the overlaps between different 
states, defined as: q ab = ^E/ij) \xtj ~ Xiji According 
to this definition, q a b = if the bosonic field has the same 
configuration in the two states, whereas q a b > otherwise. 
As in spin glasses, one can define the overlap distribution 
P( a ) = J2a,b w a w b$(q - q a ,b) where w a is the thermody- 
namic weight of the amorphous state a [6]. The overlap 
distribution is apparently continuous. P(q), averaged over 
16 different realizations of the graph is plotted in Fig. 2, 
at zero temperature and for t = 0.23. 
The transition from the uniform phase to the valence bond 
glass is continuous: the free energy of the two phases coin- 
cide within our numerical accuracy on the line t s (T) where 
the liquid phase becomes unstable. Close to the transi- 
tion point, the distribution of the Xij is peaked around 
the value x which characterizes the uniform phase, and 
it gets broader and broader as the temperature and/or 
the hopping amplitude are decreased. This transition 
shares many common features with the transition from 
the paramagnetic phase to the spin glass phase observed 
in mean field (classical) spin glasses such as, for instance, 
the Sherrington-Kirkpatrick model [6] : in both cases, one 
finds a continuous transition with a continuous distribu- 
tion of the overlaps. As a consequence it is natural to in- 
vestigate whether the VBG phase is marginally stable as 
the spin glass phase [6]. This means that the VBG phase 
is critical not only at the transition but in the whole region 
of the phase diagram where it exists. In order to do that 
we study whether the spatial correlations among valence 
bonds on different links of the lattice (xij (w n ) Xki(^n))c 
are long-ranged (as previously we focus on tu n = which 
is expected to give the main contribution). The inverse of 
the free energy Hessian matrix gives directly the dimer- 
dimcr correlations. Instead of inverting this matrix, we 
follow a less computational demanding route using a kind 
of fluctuation-dissipation relation. The idea is to measure 
the response of the system, more precisely of the value of 
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Fig. 3: Valence bond non-linear susceptibility, xvbg, as a func- 
tion of T — T a at fixed t = 0.1 (left panel) and as a function 
of t — t, at zero temperature (right panel). Xvbg diverges as 
a power law as the transition to the VBG is approached. In 
both cases the exponent is compatible with 7 ~ 1. The data are 
averaged over 8 different realization of graphs with M = 512 
sites. 

Xij , to an external perturbation and relate it to the VBG 
susceptibility. The relevant perturbation for the present 
case is a local increase of the hopping amplitude on a given 
link of the graph, t — > t+5tki- Simple integrations by parts 
in the functional integral defining the partition function, 
Eq. (3), allow one to establish the following identity: 

XVBG = ^ £ {x%Xll)l = ^ E HIT) • 

(15) 

where x% is a short-hand notation for Xijl^n = 0) and 
the subscript c denotes the connected correlation func- 
tion. We measured the response functions in the RHS of 
eq. (15). We found, as shown in Fig. 3, that the valence 
bond glass non- linear susceptibility, xvbg, diverges as a 
power law both at fixed t as the temperature is decreased 
(X2 ~ (T-T s )-i), and at fixed T (included T = 0) as the 
hopping is decreased (x2 ~ (t — £ s )~ 7 )• The exponents 
have the mean field value 7 ~ 7' ~ 1. Furthermore we 
find that xvbg is infinite (meaning of the order of, and 
scaling as, TV) in all the VBG phase, hence, confirming 
the marginality of the VBG phase. 

Differently from the transition from the liquid phase to the 
VBG, the transition from the dimer phase to the glassy 
one is discontinuous. It takes place at t c (T), where the 
free energies of the two phases coincide (at T = we have 
that t c ~ 0.175). The dimer phase becomes unstable only 
for larger values oft. Furthermore the non-linear suscepti- 
bility, Xvbg: stays finite approaching VBG from the dimer 
phase as it is expected for a first order transition. 

In summary the Valence Bond Glass phase is character- 
ized by an amorphous arrangement of dimers and absence 
of magnetic ordering. It has long-range critical dimer- 
dimcr correlations in the whole VBG phase (not only at 
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the transition). It has no gap in the electronic and spin 
density of states, although we observe a pseudo-gap. As 
a consequence it is related to, but quite different from, 
valence bond crystal and spin glass phases. We expect 
the VBG phase to be generically one of the possible low 
temperature phases arising from the interplay of strong 
quantum fluctuations and frustration. In the future it 
would be important to go beyond the simplifying frame- 
work we focused on. The role of 1 /N corrections should be 
elucidated. Furthermore, it would be interesting to study 
different models, different (and more realistic) lattices and 
add some kind of local quenched disorder. The large N 
approximation and the type of lattice we chose favor the 
glassy phase. In reality we expect that VBG will emerge 
as a true thermodynamic phase only in presence of some 
kind of quenched disorder (not much if there is already 
geometrical frustration). In this case the VBG phase will 
be in competition with the spin glass phase which in our 
treatment is excluded from the beginning because of the 
type of large N limit we used. Another interesting route 
to follow is to study the effect of doping and the result- 
ing properties of the VBG phase. This could be relevant 
for undcrdoped high T c superconducting materials. Al- 
though the VBG phase may not be a true thermodynamic 
stable phase it could nevertheless capture some kind of 
metastablc slow and glassy dynamics which seems indeed 
to be present [20] . From a more fundamental and technical 
point of view obtaining a complete solution of our model 
(analytically or by numerical simulations) would be im- 
portant to determine whether, as our results suggest, the 
VBG phase is completely analogous to the mean-field spin 
glass phase [6]. Finally it is worth studying the effect of 
magneto-elastic couplings. Because of the marginal sta- 
bility of the VBG phase they could play a very important 
role. We expect as experimental signature of the valence 
bond glass phase spatially heterogeneous NMR signals. 
Furthermore, approaching the (continuous) transition to- 
ward the VBG phase, the VBG susceptibility diverges and 
this could lead to anomalous (even divergent) non-linear 
pressure responses. Finally, we point out that preliminary 
results on modified random lattices (e.g., random regu- 
lar graphs where each site is replaced by square plaque- 
ttes) show that also glassy flux phases [10] might appear. 
These are characterized by amorphous circulating micro- 
currents. 

* * * 

We thank J. -P. Bouchaud, C. Chamon, A. Lefevre, M. 
Mezard, G. Misguich and E. Vincent for many useful and 
helpful discussions. 

Appendix. — Here we describe in detail the Monte 
Carlo annealing simulations we used. We pick up a link 
(ij) at random out of the zAf/2 total links and try to 
change either the real or the imaginary part of Xij by a 
random amount 5 G (—S max , S max ) with probability 1/2 



respectively 7 . Then we compute the new free energy, ac- 
cording to Eq. (8). Since F(\) contains a non-local term, 
at each step we have to diagonalize the matrix tC + x 
and compute all its eigenvalues, which takes a computa- 
tional time proportional to J\f 2 . The move is accepted 
with probability p = min {1, cxp[— AF/T aux ]}. The value 
of 5 r nax is self-adapted during the simulation in such a 
way that the average acceptance rate of the moves is 0.3. 
We have checked that several different values of the cho- 
sen acceptance rate lead to the same results. The aux- 
iliary temperature is decreased at constant rate down to 
very low temperature, starting from T aux = 0.5. Most of 
the results presented here have been obtained with a rate 
T = T aux /T aux = 1.3 • 10~ 3 (where each MC step consists 
of zM total attempts) . We have verified that slower cool- 
ing rates down to T ~ 5 • 10 -5 do not change the results. 
Some MC steps are finally performed at T aux = 0. 
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7 Equivalently, at each step one can also attempt to change cither 
the norm of the valence bond, |xij| 2 i by an amount 8, or its angle in 
the complex plane 9 = tan - 1 [Q(xij)/5R(xij)], by randomly choosing 
a new angle 8 in the interval (0, 27r). 



